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Abstract 

In this paper, we present a parallel higher-order boundary integral method to solve the linear 
Poisson-Boltzmann (PB) equation. In our method, a well-posed boundary integral formulation is 
used to ensure the fast convergence of Krylov subspace linear solver such as GMRES. The molecular 
surfaces are first discretized with flat triangles and then converted to curved triangles with the 
assistance of normal information at vertices. To maintain the desired accuracy, four-point Gauss- 
Radau quadratures are used on regular triangles and sixteen-point Gauss-Legendre quadratures 
together with regularization transformations are applied on singular triangles. To speed up our 
method, we take advantage of the embarrassingly parallel feature of boundary integral formulation, 
and parallelize the schemes with the message passing interface (MPI) implementation. Numerical 
tests show significantly improved accuracy and convergence of the proposed higher-order boundary 
integral Poisson-Boltzmann (HOBI-PB) solver compared with boundary integral PB solver using 
often-seen centroid collocation on fiat triangles. The higher-order accuracy results achieved by 
present method are important to sensitive solvation analysis of biomolecules, particularly when 
accurate electrostatic surface potentials are critical in the molecular simulation. In addition, 
the higher-order boundary integral schemes presented here and their associated parallelization 
potentially can be applied to solving boundary integral equations in a general sense. 
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1 Introduction 



Molecular modeling is a rising interdisciplinary approach on the study of structure, function and 
dynamics of molecules with biological significance [IJ. Among interactions in molecular modeling, 
electrostatics are critical due to their ubiquitous existence. Meanwhile, electrostatics are expensive 
to compute as they are long-range and pairwise interactions. The Poisson-Boltzmann (PB) model 
is an effective approach to resolve the electrostatics including energy, potential and force of solvated 
biomolecules [2]. As an implicit solvent model, the PB model considers solvent effects with a mean 
field approximation, and models the mobile ions with the Boltzmann distribution. These implicit 
treatments of solvent surroundings make the PB model computationally more efficient compared with 
explicit solvent models, in which atomic details of solvent molecules and electrolytes are described. 
Recently, experimentalists also showed interests in using the PB model to provide references for 
newly released structures of biomolecules e.g. the neurotransmitter receptor Acetylcholine |3j and the 
Circadian clock complex CL0CK:BMAL1 g]. 



(a) (b) 




Figure 1: (a) Poisson-Boltzmann (PB) model: domains fii (molecule) and (solvent) with different dielectric 
constants £i and £2 respectively; (b) the molecular surface is formed by the trace of solvent probe in contact 
with the solute (molecule). 



In the PB model, as illustrated in Fig. [TJ^a), the computational domain is divided into the 
solute (molecule) domain J7i and the solvent domain by a closed molecular surface F such that 

= r^i ur22 UF. As shown in Fig. [Hb), the molecular surface F is formed by the traces of a spherical 
solvent probe rolling in contact with the van del Walls balls of the solute atoms [5j, i6j . The solvated 
molecule, which is located in domain Vti with dielectric constant ei, is represented by a set of Nc 
point charges carrying Qi charge in the units of Cc at positions Xj,i = 1, ...,Nc. The exterior domain 
contains the solvent with dielectric constant 82, as well as mobile ions. For x = {x,y,z), the PB 
equation for the electrostatic potentials in each domain is derived from the Boltzmann distribution 



2 



and Gauss' law and has the form 

AT, 

V-(ei(x)V(/<i(x)) = -^gi<5(x-x,) in f^i, (1) 

j=i 

V • (£2 (x) V02 (x) ) - sinh (t>2 (x) = in , (2) 
A. ( \ A. I \ d(t)i{^) 5(/)2(x) 

0i(x =(?i2 X , ei— =£2^ on r, (3 

Oh' Of 

hm 02(x)=O, (4) 

|x|— ^-CXD 

where and 4)2 are the electrostatic potentials in each domain, qi = ecQi/ksT^i = l,...,Nc, Cc the 
electron charge, ks the Boltzmann's constant, T the absolute temperature, 6 the Dirac delta function, 
K the Debye-Hiickel parameter, and i' the unit outward normal on the interface F. We assume weak 
ionic strength in this context therefore the non-linear sinh function term can be approximated by its 
linearized term, resulting in the linear PB equation with sinh(/)2(x) term replaced by 02 (x) in Eq. ([2]). 

The linear PB model is an elliptic equation defined on multiple domains with discontinuous coef- 
ficients across the domain interfaces. The PB equation has an analytical solution only for the simple 
geometries such as spheres [7] or rods [8j- For molecules with complex geometries, the PB equation 
can only be solved numerically, which is challenging due to the following numerical difficulties. 

(1) The solutions to the PB equation, physically the electrostatic potentials, are not smooth across 
the interface as the continuities of both the potentials and the fluxes in Eq. ([3]) across the interface F 
are required to be satisfied. 

(2) The complex geometry of the interface needs to be captured to maintain the accuracy of the 
potentials particularly on or near the interface. 

(3) The partial charges carried by the individual atoms of the solute, modeled by the weighted sum- 
mation of the Dirac delta functions, is hard to accurately discretize. 

(4) The PB equation is defined on the entire domain subject to boundary condition that the 
potentials approach zero at infinity thus a cutoff for 3D mesh-based methods is inevitable. 

The wide application of the PB model as well as its associated numerical difficulties attracted 
attention from various computational science communities ranging from biophysics, biochemistry, 
mathematics, computer science, mechanical engineering as well as electrical engineering. Many nu- 
merical PB solvers were developed and they can be roughly but not completely divided into two 
categories: The 3D mesh-based finite difference/finite element methods [9l [lOl [TTl [121 US E]; and 
the boundary integral methods [l5l[ll[l7l[Il[T8l[20l[2ll[22l[23l[24]. All these methods have their 
own advantages and disadvantages. For example, the PB solvers embedded in molecular modeling 
packages such as Dephi [9], CHARMM [10], AMBER [TT], APBS [T2] use standard seven-point finite 
difference with approximated approaches to bypass the numerical difficulties (l)-(4). Although ar- 
guably these solvers have reduced accuracy, the efficient, robust and user- friendly features of these PB 
solvers brought their popularities among the bio-oriented community. An often overlooked drawback 
of these solvers is that they provide acceptable accuracy of resolved electrostatics potentials away 
from the interface but are unable to provide accurate solutions near or on the interfaces. Applica- 
tions such as molecular simulation related to ion channels [25] , cell membranes [26], and chromatin 
packing [27] require accurate electrostatic potentials and fields near or on the interface, thus call for 
developing higher-order methods to solve the PB equation. Three-dimensional mesh-based Inter- 
face methods such as Immersed Interface Methods (IIM) [13] and Matched Interface and Boundary 
Poisson-Boltzmann (MIBPB) solver [Ij] can significantly improve the accuracy by rigorously treating 
the numerical difficulties (l)-(3). However, these methods need to cope with numerical difficulty (4) 
and the complexities of the algorithms often reduce the efficiency. 

Compared with 3D mesh-based methods, the boundary integral methods have many advantages. 
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(1) The solution is characterized solely in terms of surface distributions so there are fewer unknowns 
in comparison to methods that discretize the entire domain. 

(2) The far-field boundary condition in Eq. (HD is exactly imposed. 

(3) The surface geometry of T can be represented to high precision using appropriate boundary ele- 
ments. 

(4) The electrostatic potential at charge sites are accurately determined using exact analytical ex- 
pressions. 

(5) The continuity conditions in Eq. ([3]) are explicitly enforced. 

Due to these advantages, boundary integral PB solvers have gained increased attention and we briefly 
review studies relating to the present work. In 1988, Zauhar et al. [28j introduced the boundary 
integral formulation by solving the Poisson equation for the induced surface charges. In 1990, Yoon 
et al. [29] formulated an ill-posed integral formulation of the PB equation. In 1991, Juffer et al. |15) 
reformulated the previous work to obtain a well-posed formulation, which were applied by most of the 
boundary integral PB solvers after that. The boundary integral methods can analytically circumvent 
the numerical difficulties (l)-(4), and accelerate the solver with fast algorithms such as fast multipole 
method (FMM) [HI |T71 [19] and treecode [20J . These boundary integral PB solvers mostly applied 
centroid collocation methods on fiat triangle and benchmark tests on spherical cavities with available 
analytical solutions show 0.5th order accuracy relative to number of elements [171 120] . which left 
spaces for the more challenging problem of developing higher-order boundary integral PB solver. 

In this paper, we present a more accurate boundary integral PB solver on curved triangles with 
higher-order quadratures and regularization of singularities. The rest of the paper is organized as 
follows. In section 2, we provide our algorithms including the well-posed boundary integral formulation 
and the higher-order numerical schemes, followed by the MPI parallelization. In section 3, we provide 
the numerical results, first on the spherical cavities with centered and eccentric partial changes whose 
analytical solutions are available and then on a protein (PDB: lajj) for the electrostatics solvation 
energy computation and the parallel efficiency. This paper ends with a section of concluding remarks. 

2 Methods 

We will use the well-posed boundary integral formulation from Juffer's work [15] together with high 
order quadrature [.28j and singularities regularization by a transformation [301 El] ■ We modify and 
improve these methods as needed and we will explain the details in this section. One fact affecting 
the accuracy of boundary integral methods is the discretization of the surface. For the sphere, 
we use a non- uniformed triangular surface from MSMS [32] with radial projection to correct the 
truncated output, or a uniform icosahedral triangulation [33j. For biomolecules, we only use MSMS 
to generate the fiat triangles. All fiat triangles are then converted to curved triangles by applying 
the schemes introduced as follows. Through the paper, we call our higher-order boundary integral 
Poisson-Boltzmann solver as HOBI-PB solver. 

2.1 Well-posed integral formulation 

The differential PB equation in Eqs. ([1]) and ([2]) can be converted to boundary integral equation. By 
applying the fundamental solution of Poisson equation ([T]) , Go, in r^i and the fundamental solution of 
PB equation ([2]), G K, in r22, together with Green's second theorem, and cancel the normal derivative 



4 



terms with interface jump conditions in Eq. ([3]), the coupled integral equations can be derived as [29j : 



01 (x) 



(y) 5Go(x,y) 



L(y) 



dSy + '^qkGo{x,yk), x G J^i, 



k=l 



-G«(x,y)^j— + ^- (^2(y) 



X € 



2- 



(5) 
(6) 



where Go(x,y) and G«;(x,y) are the Coulomb and screened Coulomb potentials, 



Colxjy) 



47r|x — y| 



G«(x,y) 



g-K|x-y| 

47r|x — y| 



(7) 



However, straightforward discretization of Eqs. ([5]) and ([6]) yields a linear system which becomes 
ill-conditioned as the number of boundary elements increases [34] . Juffer et al. derived a well-posed 
boundary integral formulation by going through the differentiation of the single-layer and double-layer 
potentials [15]. The desired forms are: 



+ 01 (x) 



1 + 



1\ 501 (x) 



e ) dv^ 
with the notation 



Ki(x,y)^^ + K2(x,y)0i(y) 
K3(x,y)^^ + K4(x,y)0i(y) 



dSy + s'i(x), X G r, 

(i5y -I- 52(x), X G r, 



(8) 
(9) 



i^i(x,y)=Go(x,y)-G,(x,y), 

aGo(x,y) iaG,(x,y) 



^^2(x,y) = e 



9G«(x,y) 9Go(x,y) 



^3(x,y) 



e dh 



dVy dUy ' 

52G«(x,y) a2Go(x,y) 



du^x.di'y 



du^dvy 



'S'i(x) = ^gfcGo(x,yfc), 



k=l 



(10) 



fill 



k=l 



and e = £i/£2- Note this is the well-posed Fredholm second kind of integral equation which is also 
our choice in this paper. 

In order to numerically solve the coupled equations ([8]) and Q, we need to discretize the molec- 
ular surface T with high quality elements and implement the numerical integral with higher-order 
quadrature. We also need to treat the occurred singularities or near-singularities when x and y are 
equal or nearly equal in Kernels Ki ... 4. These details are described in the following subsections. 



2.2 Curved triangles and higher-order quadratures 

Given the positions and radii of all atoms of a molecule, a triangulation program e.g. MSMS [32J can 
generate a discretized surface with a set of Nf flat triangles, Ny nodes (vertices) and corresponding 
normal directions. To achieve higher-order accuracy, we apply the schemes described in [28] to convert 
these flat triangles to curved triangles under the new background of solving the well-posed integral 
PB equation in Eqs. ([8]) and ([9]). We also modify and improve these schemes to treat the singularities 
of kernels /^i,...,4 in Eq. (jlOp . To keep this paper in an integrated form, we start from restating the 
schemes in [28] to produce the curved triangles and higher-order quadratures for discretizing Eqs. ([8]) 
and dHD. 
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Figure 3: Coordinates Transformation 



We first replace the straight element edges as shown in Fig. [2ja) with curved arcs in Fig. [2]^b) 
|28j . Let x(i) be the arc between two nodes, parameterized by the dimensionless variable t, 



x(t) 



Co + Cit + C2t^ + Cst^, 



(12) 



where co,ci, 02,03 are vector constants (12 unknowns), which will be determined by a pair of con- 
nected nodes and associated unit normals (12 conditions). At any point on the curve, the normal can 



be found by n(t) = sgn(t)- 



K{t) 



Mm 

along the curve and K{t) is the curvature given as |28j : 



, where sgn(t) = ±1 is chosen to keep a constant orientation of n(t) 



dx/dt 



dt"^ 



IP 



|dx/(it|2 dt 



(13) 



As shown in Fig. El^c), we can use parameter u € [0, 1] for the two curves starting from Xi and 
ending at X2 and X^. Then for any given n, two points on the curves X1X2 and XiX^ are specified, 
say Yi and Y2 with normal directions. By following the same procedure, we could find a curve 
connecting Yi and Y2, using v G [0, 1] as the parameter. The pair of parameters {u,v) will therefore 
correspond to a point on the curve element. 

In order to conveniently use quadrature rules, the integral is conducted in a unit right triangle. 
A mapping between the parameters r, s G [0, 1] on the unit right triangle in Fig. El^a) and parameters 
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u,v on a unit square in Fig. E^b) is constructed with the following transformation |28) . 

J u = {r + s),v = s/{r + s) if r + s 7^ , , 

\ u = 0,v = ifr = s = ^ ^ 

Based on this, we can establish the mapping from point (r, s) on the unit right triangle in Fig. El^a) to 
point x^^\u{r, s),v{r, s)) = x.^^\r, s) on the jth curved elements in Fig. ^c) by the following steps: 

(1) Given (r, s), find u,v through Eq. ST^ ■ 

(2) Plug u as the parameter into the curve functions of X1X2 and X1X3 to locate Yi and Y2 as shown 
in Fig. EJc). 

(3) Find the curve function of Y1Y2 and then plug v as the parameter to finally get x^-^^ on the curved 
elements. 

However, this is not an analytical function for efficient computations but a constructive procedure. 
As a remedy, a high-order 10-point interpolation scheme is used [28]. The brief idea is illustrated in 
Fig. ([3]) on the jth element of the triangulation. 

(1) Pick 10 specified points (r^, Sk) first in the right unit triangle for k = 1,2, 10. 

(2) Find u{rk, Sk),v{rk, Sk) according to Eq. (fH|) . 

(3) Find x-^(n(rfc, s^), u(rfc, s^)) and their normal directions on the curved element, which is param- 
eterized by {u,v). Note points 1,4,9 are already given with the flat triangles. With these three 
nodes, three trajectories can be found and used to find the positions and normal directions of points 
2,3,5,6,7,8. Finally with points 5 and 6, the dashed trajectory can be formed to find point 10 and its 
normal direction. 

(4) As required by the quadrature rules, interpolate any point on the curved element by the 
expression 

10 

^(3\r,s)=Y,Nk{r,s)^^P (15) 
k=l 

where Nk{r, s) are Lagrangian interpolation polynomial for a 10-point element. See table 1 of [28j for 
the expression of Nk{r,s). Note Eq. (fTSll can also be used to find the partial derivative of x^-^^ with 
respect to r and s, which is required for computing the Jacobian for transformation and the normal 
direction at x^-'^. 

Suppose now we will integrate function /(x) on the jth curved element. The quadrature rules 
give the position of a set of points on the unit right triangle e.g. (rm, Sm) and quadrature weights Wm 
for m = 1, 2, .., N, where is the number of quadrature points, and the integral can be evaluated as: 

Note the term ^^^^^'^""^ x ^^^'g^'^™'^ gives the normal direction of the point x^^\r, s) with parameters 
(r, s), and we will use this information to supply the required normal direction in Eqs. ([SD and Q- 

In this paper, we use N = Nqr points (practically we choose Nqr = 4) in the Gauss-Radau 
quadrature [35]. For i = 1,2, ...,Ny, the ith and the (i + A'^)th element of the discretized matrix- 
vector product Au are given as 



N 



9x(rv^^mJ 9x(r. 
i))\ ~ ^ 



771=1 



dr 



ds 



IWrr 



(16) 



, K Ngr 3 

j — 1 m—l n—l 



Ngr 3 



j — 1 m — l n — l 



—Q^ + i^4(xi,y^-,„)<;!>i(y^„) 



(17) 
(18) 
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In Eqs. ()17p and (jlSp . the superscript r of {Au} stands for regular triangle, Wj^rn,n contains the weights 
and coefficients associated with the quadrature, the transformation Jacobian, and the interpolation 
coefficients in Eq. (jl5p . Note that is the number of regular triangles associated with the ith vertex. 
In these equations, y^^^ are the same set of nodes as Xj, but Yj^rn quadrature points on the jth 

curved element, which are mapped from the predetermined points on the unit right triangle. This 
mismatch brings difficulty to apply Fast Multipole Method or treecode to accelerate the higher-order 
scheme. Studies about this issue will be proceeded in our future work. 



2.3 Treatment of singularities 

From Eqs. (]17p and (jlSp . we can see that when node i is one of the vertices of the jth curved triangular 
element, singularity or near singularity occurs at evaluating kernels i^i,...,4. In other words, Xj is 
equal to or nearly equal to yj,m- It can be shown that the singularities in kernels -fCi,...,4 are in the 
order of 0( i — r) 1 151. To treat these singularities, we use the tensor-product of Gauss-Legendre 
quadrature on a unit square, together with a transformation [301 131j . In this case the mapping from 
a unit square (0 < x, ?/ < 1) to a unit right triangle (0 < r, s < 1 and r + s < 1), and to a curved 
triangle is constructed. The mapping that r = (1 — y)x and s = for < x, y < is used to 
remove the singularities in kernels i^i,...,4 when they appear at (r, s) = (0, 0) which indicates x = 0. 
To understand this, it can be seen that the Jacobian for the transformation from (r, s) to (x, y) is x, 
thus it can remove the 0( 1 _\ 1 ) type of singularities. 

I^z yj,m 

In the treatment of singularities, if the number of quadrature points used in each direction of the 
Gauass-Lagendre quadrature is Nql, the ith and {i + N^)th. element of the discretized matrix- vector 
product Au will in addition contain the following singular component 



{AuK = ^{l+e)M^^)^Yl E E^^^-.-." 

j—1 m—l n—l 



' +^2(xi,y. )(;/)i(y, ) 



1 (NoLr 3 



'Wi. 



j—1 m—l n—l 



J«'3(x»,y. 



) +-^4(x.,y^.^)0i(y^.„) 



(19) 
(20) 



In Eqs. (fTop and (pO]) . the superscript s of {Au} stands for singular triangle, and Wj^rn,n contains the 
weights and coefficients associated with the quadrature, the transformation Jacobians (additionally 
contains the Jacobian of the r = (1 — y)x and s = yx mapping), and the interpolation coefficients 
in Eq. (jlSp . Note index j has the range up to Nf , which stands for the number of singular triangles 
associated with the ith vertex. Simulation shows this value is various for different vertices and it 
can be as big as 15. In this paper, we practically choose Nql = 4 points in each direction to ensure 
desired accuracy [35] . 



2.4 Low-order scheme 

We also briefly introduce the low-order boundary integral Poisson-Boltzmann (LOBI-PB) solver, 
which is used for comparison in this paper. In this low-order scheme, the flat triangle and centroid 
collocation are used, i.e. the quadrature point is located at the center of each triangle. This scheme 
also assumes that the potential and its normal derivative, as well as the kernel functions are uniform 
on each triangle. When singularity in kernels occurs, the contribution of this triangle in the integral 
is then simply removed. This scheme is in fact widely used in latest boundary integral methods in 
solving PB equations to provide convenience on incorporating fast algorithms such as EMM [34J and 
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treecode [20]. For i = 1,2, Nf, the ith and the {i + A'^y)th element of the discretized matrix-vector 
product Au are given as 



{An}i = -(l + £)0i(x,)- 

Nf 



OV-v 



1 



1 



L(xi) 



(21) 



(22) 



It worths mentioning that the unknowns of LOBI-PB is at the centroid of the triangle in the number 
of Nf while the unknowns of HOBI-PB is at the vertices of the triangle in the number of N^. 



2.5 Electrostatic Solvation Energy Formulation 

The electrostatic solvation energy is computed by 

A^c , N, 



Esol = ^ E 9fc'/'reac(xfc) = ^ ^ ^'^ / 
k=l k=l 



Ki(x,,y)^^+K2(xfc,y),/.i(y) 



dSy, 



(23) 



where (j)rea.c{^k) = (pii^k) — 5'i(xfe), whose formulation is the integral part of Eq. (p3]) . is the reaction 
potential at the A;th solute atom. The electrostatic solvation energy, which can be regarded as the 
atomistic charge weighted average of the reaction potential </)reac5 can effectively characterize the 
accuracy of a PB solver. 



2.6 MPI Based Parallel Implementation 

The HIBO-PB solver uses the boundary integral formulation, which can be conveniently parallelized. 
The majority of the CPU time is taken by the following routines. 

(1) Compute the source term in Eqs. ^ and ([9]). 

(2) Convert flat triangles to curved triangles. 

(3) Compute and store the Gauss-Radau quadrature information for each curved triangle. 

(4) Compute and store the Gauss-Legendre quadrature information for singular triangles associated 
with each vertex. 

(5) Perform matrix- vector product as in Eqs. ()17p . (Iisp . (jl9p . and ()20p on each GMRES iteration. 

(6) Compute electrostatic solvation energy with Gauss-Radau quadrature. 

Among all these routines, routine (5) is the most expensive one and it will be repeatedly computed 
on each GMRES iteration. We parallelize all these routines to maximize the parallel efficiency. Table 
[1] provides the Pseudocode for MPI based parallel implementation. 



3 Results 

In this section, we present numerical results. We first solve the PB equation on spherical cavities 
with a centered charge and with an eccentric change at different locations. The analytical solutions in 
terms of a closed form (for centered charge) and in terms of spherical harmonics (for multiple eccentric 
charges) are available for these tests [7j. To demonstrate the higher accuracy obtained by higher- 
order boundary integral Poission-Boltzmann (HOBI-PB) solver, we compare the numerical results 
with APBS[12j, MIBPB[14], and the lower-order boundary integral Poisson-Boltzmann (LOBI-PB) 
solver. APBS uses straight-forward finite difference scheme. MIBPB is a 2nd order interface method 
repeatedly using local interpolation to capture interface jump conditions [Ml [37j and applying a 
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Table 1: Pseudocode for parallel HOBI-PB solver using replicated data algorithm. 



1 on main processor 

2 read protein data 

3 call MSMS to generate triangulation 

4 copy protein data and triangulation to all other processors 

5 on each processor 

6 locally compute sources terms for each assigned vertex 

7 locally convert flat triangles to curved triangles 

8 for each assigned triangle, locally compute and store quadrature information 

9 for each assigned vertex, locally store quadrature of associated singular triangles 

10 copy result to all other processors 

11 set initial guess for GMRES iteration 

12 compute assigned segment of matrix- vector 

13 copy result to all other processors 

14 on main processor 

15 test for GMRES convergence 

16 if no, go to step 12 for next iteration 

17 if yes, go to step 17 

18 on each processor 

19 compute assigned segment of electrostatic solvation energy 

20 copy result to main processor 

21 on main processor 

22 add segments of electrostatic solvation energy and output result 



Dirichlet-to- Neumann mapping to transform the singular charges to interface jump conditions |38j . 
LOBI-PB discretizes the integral equations with flat triangles and performs the numerical integral 
with centroid collocation as explained in the previous section. 

We then solve the PB equation on a test protein (PDB: lajj), which is a lipprotein receptor 
with 37 residues and 519 atoms. We report the electrostatic solvation energy results for this protein 
computed from both HOBI-PB and LOBI-PB to demonstrate the improved accuracy achieved by the 
higher-order integral schemes. All algorithms are written in Fortran 90/95 and compiled with GNU 
Fortran with flag -03. The serial simulations are performed with a single CPU (Intel(R) Xeon(R) 
CPU E5440 @ 2.83GHz with 2G Memory) on an 8-core workstation. MPI parallel simulations are 
conducted on the DMC cluster (Intel(R) Xeon(R) CPU E5520 @ 2.27GHz with 1.5G memory for 
each core) at Alabama Supercomputing Center. Before seeing the numerical results, we define order 
and errors. 

3.1 Order and Errors 

In this paper, we report the relative errors defined as 

max I (/.""'"(xi )-(/''='"' (xi) I 

^ i = l,...,A'' (24:) 

max \<h'^^°-(xi)\ 
i=l,...,N 

where N is the number of unknowns. Note is the number of vertices for HOBI-PB, the number of 
triangular elements for LOBI-PB, and the number of close-to-surface mesh points (irregular points) 
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Table 2: Accuracy comparison of different PB solvers on a spherical cavity (radius=2A , (7(0,0,0) — ICc, 
e = 80, K — Q)\ h the mesh size; d the number of vertices per A^. 





APBS 


MIBPB 




HOBI-PB 


LOBI-PB 


h 


Esol 


60 


ord. 


Esol 


60 


ord. 


d 


Esol 


60 


ord. 


Esol 


60 


ord. 


1 


-83.44 


1.94e+0 




-81.95 


1.24e-2 




5 


-81.98 


1.45e-4 




-83.41 


4.07e-4 




0.5 


-85.85 


1.31e-|-0 


0.6 


-81.98 


1.91e-3 


2.7 


10 


-81.98 


5.26e-5 


1.5 


-83.13 


2.55e-4 


0.7 


0.2 


-82.58 


5.76e-l 


0.9 


-81.98 


3.87e-4 


1.7 


20 


-81.98 


1.52e-5 


1.8 


-82.82 


1.82e-4 


0.5 


0.1 


-82.27 


2.94e-l 


1.0 


-81.98 


1.07G-4 


1.9 


40 


-81.98 


6.13C-6 


1.3 


-82.60 


1.27e-4 


0.5 


0.05 


-82.03 


1.49e-l 


1.0 


-81.98 


2.31G-5 


2.2 


80 


-81.98 


1.85C-6 


1.7 


-82.43 


8.58e-5 


0.6 



for APBS and MIBPB. The notation 0"-"™- represents numerically solved surface potentials and (j)'^^"' 
denotes the analytical solutions obtained by Kirkwood's spherical harmonic expansion \J\. The dis- 
cretization of APBS and MIBPB are on the Cartesian grid with mesh size h. The discretization of 
HOBI-PB and LOBI-PB are on the molecular surface with density d, number of vertices per A^. 
The numerical order of accuracy is computed with 

coarse_error 

order = log coarse_mesh (25) 

fine.mesh fine_error 

following the convention of numerical analysis, where "mesh" refers to h for finite difference methods 
or density d for boundary integral methods, both at coarse and fine levels. 



3.2 On a Spherical Cavity with one Centered Charge 

We first solve the linear Poisson-Boltzmann equation on a spherical cavity with radius 2A and a 
centered change lec submerged in water with zero ionic strength. We report the electrostatic solvation 
energy Egoi and surface potential errors e^f, computed with all above-mentioned methods in Table [2j 
The results of APBS and MIBPB are from reference [38]. 

From Table [21 we can see that APBS provides acceptable value in electrostatic solvation energy 
compared with the true value 81.98 kcal/mol since only the potential at the center of the spherical 
cavity is required to compute the electrostatic solvation energy. However, the surface potentials com- 
puted by APBS method show large errors. This is due to the approximation on interface conditions 
and singular charges of standard finite difference methods. 

MIBPB uses a more sophisticated finite difference scheme. The rigorous treatment on interface 
conditions and singular charges significantly improves accuracy. The electrostatic solvation energy 
is nearly perfect even at coarse grid and the surface potential is very accurate with solid 2nd order 
convergence pattern relative to mesh-size h. 

LOBI-PB gives the electrostatic solvation energy in the same accurate level as APBS but produces 
much more accurate surface potentials than APBS does. These surface potentials converge at the 
0.5th order relative to density d. 

HOBI-PB is obviously a more accurate method. It shows 1.5th order of convergence relative to 
density d as reflected from surface potential errors and its accuracy is even better than results from 
MIBPB. 

These results demonstrate HOBI-PB is the most accurate PB solver among its peers. According 
to our knowledge, for this classic benchmark test, there is no other PB solvers can achieve the same 
level of accuracy as HOBI-PB does. 

3.3 A Spherical Cavity with One Eccentric Charge 

We further investigate the performance of HOBI-PB and LOBI-PB on a spherical cavity with radius 
lA and a lec eccentric charge at different locations. The errors shown are in terms of surface potential 
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Cu as previously defined. The Debye-Hiickel parameter k is set to 1. The charge moves from the 
center of the sphere toward to the surface and we test the performance of the PB solvers in response to 
the change of locations. The closer the charge is to the surface, the more variant the induced charges 
on the surface appear to be. This interesting phenomenon draws attention from many researchers, 
e.g. the peak separation method by Juffer et al. |15j and the image method by Deng et al. ^39j. In 
addition to that, we also try to investigate if the quality of the triangulation will affect the accuracy. 
To this end, we report the results computed on the triangular surfaces generated by MSMS [32] and 
by a uniform icosahedral triangulation routine (ico) |33j . MSMS generates triangles in various shapes 
and sizes, while ico generates uniform and equilateral triangles with fixed numbers of triangles such 
as 20, 80, 320, etc. The results are plotted in Fig. HI and we have the following observation and 
discussion. 

(1) The results of HOBI-PB in Fig. (b)(d) show significant improvement in accuracy compared with 
results of LOBI-PB in Fig. (a)(c). For HOBI-PB, the errors are smaller in general and the difference 
between two different meshes are larger, indicating better accuracy and higher-order convergence 

(2) LOBI-PB, due to its simplicity in algorithm, shows more consistent convergence pattern. HOBI- 
PB, however affected by the quadrature rule and singularity removal schemes, shows some irregular 
patterns. 

(3) We also observe that when the charge is located close to the interface e.g. at (0.9,0,0), for coarse 
mesh, LOBI-PB shows better accuracy than that from HOBI-PB. This is due to the fact that the 
close-to-boundary charge brings high variation of the induced charges on surface and multiple quadra- 
ture points selected in HOBI-PB amplify the variations. In LOBI-PB, there is only one quadrature 
point in each triangle at the centroid thus the scheme is less sensitive to the induced charges. 

(4) The quality of triangular surface will slightly affect the convergence of the boundary integral PB 
solvers. The results from icosahedral triangulation routines show more uniform convergence pattern 
than that from MSMS. 

(5) We further draw the error-vs-element plots for one charge located at (0,0,0) in Fig. \M,e) and 
(0.9,0,0) in Fig. S^f) for both LOBLPB (circle) and HOBI-PB (triangle) solvers on both MSMS (red, 
empty marker, solid line) and icosahedron (blue, filled marker, dashed line) triangles. By observing 
one color at a time (one kind of mesh at a time) , we can see the pattern in Fig. W^e) is uniform and 
HOBI-PB shows better accuracy (smaller y values) and faster convergence (larger slope). The pattern 
in Fig. m^f ) is tangled. By observing the slope, we can still see that HOBI-PB converges faster than 
LOBI-PB generally. However, the error of the HOBI-PB is bigger than LOBI-PB for coarse mesh due 
to the interaction between the induced charges on surface and the singular charge near the surface, as 
explained in observation and discussion (3). In practice, the partial charges are at least 1-2 A away 
from the molecular surfaces therefore the slightly deteriorated pattern observed here when charges 
are close to the surfaces will unlikely happen. 

In short, the numerical results of HOBI-PB and some other reference PB solvers on the spherical 
cavities compared with available analytical solutions quantitatively demonstrate the achieved better 
accuracy and higher-order convergence of HIBO-PB. The complexity of the higher-order schemes 
for quadrature and regularization introduces only minor instabilities but gain significantly improved 
accuracy and convergence. Next we use HOBI-PB and LOBI-PB to solve PB equation and compute 
electrostatic solvation energy on a protein. 

3.4 Computing Electrostatic Solvation Energy on Protein lajj 

The ultimate goal of HOBI-PB is to provide accurate electrostatic potentials for solvated biomolecules. 
With this solver, we solve the PB equation and compute the electrostatic solvation energy for many 
small-middle sized proteins. Here we take protein lajj with triangulation at different resolutions 
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Figure 4: relative surface potential errors (e^) on a spherical cavity with radius lA and eccentric unit charge, 
K = 1, and e = 80: (a) LOBI-PB, MSMS; (b) HOBLPB, MSMS; (c) LOBLPB, ico; (d) HOBLPB, ico; (e) 
error vs. elements ^ N, a charge located at (0,0,0); (f) error vs. elements # iV, a charge located at (0.9,0,0). 
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as an example. The coordinates and partial charges of the protein are obtained from CHARMM22 
force field [IQ]. The numerical results show that HOBI-PB solves the PB equation with significantly 
improved accuracy compared with LOBI-PB does. 

Table 3: Testing results for protein lajj. HOBI-PB results showing electrostatic solvation energy Esoi, CPU 
time, memory usage; number of GMRES iterations (it.). 



d 


# of ele. N 


Esoi (kcal/mol) 


CPU (s) 


Memory (Mbyte) 


# of it. 


1 


6027 


-1168.87 


115 


47 


10 


2 


9198 


-1152.42 


266 


70 


10 


4 


17278 


-1145.50 


867 


129 


9 


8 


32386 


-1140.60 


3114 


238 


9 


16 


66558 


-1139.23 


17790 


523 


13 


32 


132028 


-1138.38 


56039 


759 


10 


64 


270680 


-1138.49 


254198 


2116 
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The numerical results for computing electrostatic solvation energy for protein lajj are reported 
in Table [3l To produce these results, we discretize the molecular surface of protein lajj at different 
densities as seen in the first column of the table. Different densities result in different numbers of 
triangular elements which are listed in the second column of the table, characterizing the dimension 
of problem. In the third column, we report electrostatic solvation energies. We can see these values 
are very close to each other at different resolutions and they converge to a value near about —1138.49 
kcal/mol. In column 4, we report the CPU time and it increases at the order of 0{N'^). This reveals 
currently the most critical limitation of HOBI-PB as the 0{N'^) computational cost eventually will 
make the HOBI-PB prohibitively expensive. To alleviate the pain, we take advantage of the convenient 
parallelization of boundary integral formulation, and we will see the parallelization performance next. 
Memory uses are shown in column 5 and we see a 0{N) pattern, which is advantageous compared 
with the 3D mesh-based methods whose memory uses increase at 0{N ) or even 0{N'^), where N 
is the number of unknowns. The last column is the number of iterations, and these stable results 
attribute to the well-posed integral formulation in Eqs. ([8]) and ([9]). 

We next plot the solvation energies of HOBI-PB (solid circle) from Table [3] on Fig. [5] together 
with the solvation energies computed with LOBI-PB (empty circle) for the same protein at different 
resolutions. The plot shows results from LOBI-PB converge to that of HOBI-PB. By using cubic 
interpolation, we can see both methods eventually converge toward almost the same values (the 
red "x" for HOBLPB and the red "*" for LOBLPB). The advantage of HOBLPB is that it can 
achieve high accuracy even at very coarse mesh. For example, the first solid point from right, which 
corresponds to electrostatic solvation energy of —1168.87 kcal/mol at density d = 1, is only 30 
kcal/mol different from the interpolated true value at about —1139.09 kcal/mol. Similar patterns are 
observed on our other tests on different proteins. 

We finally provided the MPI based parallel performance in terms of parallel efficiency of HOBI-PB 
for protein lajj at different meshes in Table [H Due to the limitation of the resources, our results 
are generated with up to 64 CPUs. The maximum memories for each CPU in MPI implementation 
is about 1.5G therefore for protein (lajj) we solve the PB equation with maximum d = 32 requiring 
759M memory (note d = 64 requires more than 2G memory per core from Table [3|). From Table HI we 
can see the CPU time is substantially reduced when the code is run in parallel on high performance 
computers. For example, for d = 32 with 132,028 elements, the serial work requires more than half of 
a day (50457s), the parallel work with 64 CPUs produce results in about 15 minutes (860s). We see 
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Figure 5: Comparison of electrostatic solvation energy computed with HOBI-PB and LOBI-PB for protein 
lajj 



high parallel efficiency from the column. When the dimension of the problem is sufficiently large 
(e.g. N = 66,558 or TV = 132,028), the parallel efficiency with 64 CPUs is higher than 90%. We 
observed occasionally the interesting larger-than-one parallel efficiencies and those could be explained 
by the traffic fluctuation on the cluster or possibly the argument mentioned in |41j . 

Table 4: MPI parallel performance for computing electrostatic solvation energy on protein lajj; p is number 
of CPUs, Tp is the time using p CPUs, Ti/pTp is the parallel efficiency. 





N = 


132028 


N = 


66558 


N = 


32386 


p 


Tp{s) 


Ti/pTp 


Tpis) 


Ti/pTp 


Tpis) 


Ti/pTp 


1 


50457 


100.0% 


17755 


100.0% 


3060 


100.0% 


2 


24740 


102.0% 


8790 


101.0% 


1528 


100.1% 


4 


12877 


98.0% 


4381 


101.3% 


777 


98.5% 


8 


6398 


98.6% 


2190 


101.3% 


406 


94.3% 


16 


3321 


95.0% 


1090 


101.8% 


194 


98.8% 


32 


1677 


94.0% 


571 


97.2% 


106 


90.6% 


64 


860 


91.7% 


306 


90.8% 


62 


77.4% 



In summary, the HOBI-PB solver solves the PB equation accurately on both spherical cavities 
and real biomelecules. The surface potential and electrostatic solvation energy computed with the 
solver is accurate, fast-convergent and stable. The parallelization of the solver is easy to implement 
and the parallel efficiency is attractively high. 

4 Conclusion 

This paper describes the schemes of a higher-order boundary integral Poisson-Boltzmann (HOBI-PB) 
solver. This solver discretizes the molecular surface with curved triangles, and performs numeri- 
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cal integral with four-point Gauss-Radau quadratures on regular triangles and with sixteen-point 
Gauss-Legendre quadratures on singular triangles. The singularities are regularized with a coordi- 
nate transformation. The numerical tests on spherical cavities show that HOBI-PB can achieve 1.5th 
order convergence of on surface potentials relative to boundary elements and these computed surface 
potentials are very accurate even at coarse mesh. In addition, the order of convergence does not 
compromise when the electric charge is off-center or even closed to the surface. The numerical tests 
of HOBI-PB on biomolecules show much improved accuracy compared with results from the popular 
LOBI-PB solvers. The accurate surface potentials are of vital importance to molecular modeling 
that are sensitive to electrostatics near or on the molecular surface. To improve the efficiency of 
HOBI-PB, we also developed its MPI based parallel version. The numerical results demonstrate very 
encouraging parallel efficiency, e.g. above 90% when up to 64 CPUs work concurrently. 

HOBI-PB achieves higher accuracy at the price of more complex algorithms. The limitation of 
the HOBI-PB is mainly at the problem dimension it can treat, which is subject to the available 
computing resources. For example, for the accessibility to up 64 CPUs each with 1.5G memory per 
core at Alabama Supercomputing Center, the problem size HOBI-PB can handle for a reasonable 
long waiting time (<15 minutes) is about 150,000 elements. Considering the high accuracy, we can 
use fairly small density 1 < d < 5, then PB equations on proteins with hundreds of residues can be 
conveniently solved. More computing resources will bring better performance on bigger problems. 
The rapid updating of computing power will definitely make HOBI-PB more and more capable. 

There are many spaces in which HOBI-PB can be improved and extended. For example, we are 
looking for better triangulation programs for the molecular surfaces |42 1I43[[51] . The currently adopted 
MSMS only provides 3 digits accuracy in vertices, positions and normal directions (we radially project 
vertices for spheres). In addition, the adoption of fast algorithms such as FMM [35]] and treecode [l6] 
to HOBI-PB, although considerably challenging due to the complexity of HOBI-PB schemes, is under 
our consideration. A more challenging problem is the application of HOBI-PB to molecular dynamics 
|47[|48j . where the PB equation will be solved at every time sampling. Furthermore, the higher-order 
schemes applied in HOBI-PB has the potential to solve other integral equations such as the integral 
forms of Helmholtz equation [39] and Maxwell Equations [50]. For these challenges, the application of 
the quadrature rules and the treatment of singularity will be similar, however, new challenges such as 
obtaining the well-posedness of the integral formulation and applying the fast algorithms need further 
investigation. 
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